Notes on causal differencing in ADM/CADM formulations: a ID comparison 
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Causal differencing has sfiown to be one of tlie promising and successful approaches towards 
excising curvature singularities from numerical simulations of black hole spacetimes. So far it has 
only been actively implemented in the ADM and Einstein-Bianchi 3+1 formulations of the Einstein 
equations. Recently, an approach closely related to the ADM one, commonly referred to as as 
"conformal ADM" (CADM) has shown excellent results when modeling waves on flat spacetimes 
and black hole spacetimes where singularity avoiding slices are used to deal with the singularity. In 
these cases, the use of CADM has yielded longer evolutions and better outer boundary dependence 
than those obtained with the ADM one. If this success translates to the case where excision is 
^ , implemented, then the CADM formulation will likely be a prime candidate for modeling generic 

^ • black hole spacetimes. In the present work we investigate the applicability of causal differencing to 

I CADM, presenting the equations in a convenient way for such a goal and compare its application 

. with the ADM approach in ID. 



> 

^ \ I. INTRODUCTION 

o 

. One of the goals of numerical relativity that has proven to be elusive (using a 3+1 splitting of the Einstein equations) 
has been that of modeling a generic single black hole long periods of time. Present single black hole simulations in 3D 

, have not yet been shown to be generically stable. There are limited instances of stability based on the outer boundary 

' choice and placement Most simulations run just beyond a few hundreds M based on outer boundary placement 

Q \ and binary black hole simulations run for about 20 — 50A/ before the codes either crash or the entire grid is inside the 

O^' event horizon. In some cases, the reason of the crash is well understood. For instance, the use of singularity avoiding 

5^1^ [ slices lead to the presence of steep gradients which eventually can no longer be handled by the codes. A solution to 

bJ[)- this problem is to "excise" the singularity from the computational domain Unfortunately, in most cases, it is not 

^ ■ clear what the main reasons behind the crash are and consequently addressing the problem becomes cumbersome. In 

• i-H , attempting to deal with this issue there are several possible avenues to either remove or provide an understanding of 

■ the source of problems. These avenues can be divided in the following way: 
H ' 

' 1. Choice of formulation of Einstein equations; 



2. Choice of gauge; 

3. Numerical implementations. 

Avenue (1) is motivated by the difficulties encountered in achieving long term evolutions with the ADM formulation, 
which historically has been the main tool in Numerical Relativity. Several formulations exist in the literature that 
exhibit properties like hyperbolicity the equations are expressed in a flux conservative form |5| and/or try to 
separate transverse modes Avenue (2) is based on the fact that, in principle, a coordinate system could be 

chosen such that the fields vary slowly in time; hence, the simulations would be better behaved. Conditions to achieve 
such coordinates have been presented in the literature ^ . Lastly, avenue (3) highlights the need for a more profound 
understanding of the numerical implementation of the evolution equations. Algorithms specifically tailored to deal 
with the equations under study could pave the way to better behaved simulations (for instance, compare with the 
implementations that deal with the fluid equations and their 'historical evolution' from crude implementations in 
early simulations to high resolution shock capturing schemes in present state of the art codes). 

Our present work focuses primarily on avenue (1); although avenues (2) and (3) also play a role since notable 
improvements are achieved with specific gauge choices and the use of causal differencing algorithms. We compare 
results obtained from the use of black hole excision with two related forms of the standard ADM 3+1 equations; 
focusing on the use of the CADM system with excision techniques in spherical symmetry. 
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The main motivation behind the comparison with the conformal ADM is the report by many groups that robust 
implementations have been achieved in hnearized gravity, gravitational wave spacetimes, systems containing matter, 
etc However, so far, it has only been used to model black hole spacetimes using singularity avoiding slices 

As it is widely accepted, these types of slicings are useful when the desired simulation time is rather short. 
In order to model black hole spacetimes for long periods of time, singularity excision must be employed. To study 
the feasibility of excision in this formulation and to analyze its advantages and disadvantages with respect to the 
traditional ADM formulation (where excision techniques have been used for several years already we present 

a ID study and compare results obtained with both approaches. We start with a brief review of the formulations in 
section II. In section III, we rewrite the system of equations in a way convenient for causal differencing and describe 
how this techniques is implemented. In section IV we compare simulations of a Schwarzschild black hole and show 
how the ADM formulation yields longer term evolution unless the trace of the extrinsic curvature is frozen in time, in 
which case CADM yields better behaved evolutions than the ADM formulation, we also show how causal differencing 
indeed gives the expected results in terms of stability. Wc conclude in section V with a brief discussion. 



II. FORMULATION 

The standard ADM equations, in the form most commonly used in numerical relativity, are ul 



~D,Dja + a[Ri.j + KK, 



(2.1a) 



(2.1b) 



with 



d 
dt 



dt- JOr, 



(2.2) 



where Cp is the Lie derivative along the shift vector i?^ is the Ricci tensor and Di the covariant derivative 
associated with the three-dimensional metric jij . 

The conformal ADM equations ^,|^ are obtained from the ADM ones by (I) making use of a conformal decomposition 
of the three-metric as 



= e-^'^l^J with e^"^ = 7^/3 ^ det(7,,)i/3 



(2.3) 



(hence det{j) = 1); (II) decomposing the extrinsic curvature into its trace and trace-free components. The trace-free 
part of the extrinsic curvature Kij , defined by 



(2.4) 



and K = Kij is the trace of the extrinsic curvature and (III) further conformally decomposing Aij as: 



Aij — ^ ^ Aij 

In terms of these variables, Einstein equations in vacuum are 



(2.5) 



d , 
dt"^ 



= -2aA. 



I J 1 



-aK . 
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—K = -fW,Dja + a 



Ai,A'^ + -K^ 



e-'^'^ [-DiD^a + aR,j] 
+a [kA,j - 2AuA] 



TF 



(2.6a) 
(2.6b) 
(2.6c) 

(2.6d) 
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where the Hamiltonian constraint was used to ehminate the Ricci scalar in equation (2.6c). Note that with the 
conformal decomposition of the three-metric, the Ricci tensor now has two pieces, which are written as 



Rij - Rij + 



(2.7) 



The "conformal-factor" part i?f is given directly by straightforward computation of derivatives of i 



(2.8) 
(2.9) 



while the "conformal" part Rij can be computed in the standard way from the conformal three-metric ^ij . 

To this point, the equations have been written by a trivial algebraic manipulation of the ADM equations in terms of 
the new variables. The non-trivial part comes into play by introducing what Ref. Q calls the "conformal connection 
functions" : 



(2.10) 



where the last equality holds since the determinant of the conformal three-metric 7 is unity. Using the conformal 
connection functions, the Ricci tensor is written as: 



Ri 



1 



ijl 



^ kj^ li ■ 



(2.11) 



Where are to be considered independent variables whose evolution equations are obtained by a simple commutation 
of derivatives. 



dt 



(2.12) 



As proposed in Ref. 0| the divergence of A^^ is replaced with the help of the momentum constraint to obtain: 



d ~. 
dt 



(2.13) 



With this reformulation, in addition to the evolution equations for the conformal three-metric 7^ (2.6a) and the 
conformal-traceless extrinsic curvature variables Aij ( ^.6d ), there are evolution equations for the conformal factor 
(j> (2.61), the trace of the extrinsic curvature K (2.6c) and the conformal connection functions F* (2.13). 



III. CAUSAL DIFFERENCING IMPLEMENTATION 

Causal differencing, as explained in |p^po||, provides a straightforward way to integrate the evolution equations 
while preserving (and taking advantage of) the causal structure of the spacetime under consideration. In the approach 
used in the present work we follow the strategy described in [Q. First, the Lie derivative along /3 is split and terms 
containing derivatives of P are moved to the right hand side. Then, the ADM system of equations is reexpressed as 

doltj = -2a7y + 27/(j/3^) , (3.1a) 
^oK^J =D,Dja + a + KK^, 

~2K,kK'', - (4)^^^.^ ^ 2if,(,/3^) ; (3.1b) 

and the CADM system of equations then reduces to 



3 



dolij = 




do<t> = 


















-YWiDja + a 











a KA, 



2AuA\ 



2a 



QA'- 



(3.2a) 
(3.2b) 

(3.2c) 



(3. 2d) 



(3.2e) 



where do = dt — P^di. 

Finally the numerical implementation of the equations is split into two steps. First, the equations are evolved along 
the normal to the hypersurface (at constant t) n°- — df — /3*9f . In the second step, an interpolation is carried over to 
obtain values on grid coordinate locations (see Fig [^) . Note that the two systems of equations have in this form the 
same basic structure; hence, simple modifications to an ADM code with excision (like AGAVE m^) enable the 
use of already developed excision modules with the CADM equations in a straightforward mannerff 




n+l 



FIG. 1. Illustration of the causal differencing strategy. First the integration proceeds along the dashed lines to obtain 
values in the n+l level (at filled square points). Then, an interpolation is carried out to obtain values on the grid points 
(filled circles). In the graph, as an example, a second order interpolation (indicated with arrows) provides values on the i-th 
grid point. 

A. Numerical Implementation 

In the numerical implementation of the CADM it is convenient to introduce an intermediate variable F such that 
(j) = l/41n(F) and evolve F instead of 0. This choice avoids unnecessary handling of exponential and logarithmic 
functions thus preventing loss of accuracy; hence, the equation for F is 



doF = -^aFK+^F(3^^, 



(3.3) 



A second order finite difference code has been written to implement both the ADM and CADM formulations. The 
same causal differencing algorithms have been applied to both formulations. These algorithms involve at their core, an 
interpolation for every grid point in the computational domain. Near an excision boundary the choice of interpolation 
order and stencil choice becomes important. We allow a choice of second, third and fourth order interpolations in 
order to study possible practical approaches |^. 



^ These modifications are in place in the AGAVE code and currently being tested in 3D 
■^This code is publicly available and can be requested from the authors. 
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IV. APPLICATIONS 



To compare evolutions with the above formulations we pick as a particular example the Schwarzschild spacetime 
(and linear perturbations of it). In order to implement excision, a slicing must be chosen such that surfaces of 
constant time "penetrate" the horizon. The ingoing Eddington Finkelstein coordinates define hypersurfaces satisfying 
this condition; in terms of them, the line element reads 

ds^ = - il - —jdt^ + —dtdr + f 1 + — ] + r^dn^ , (4.1) 
where dfP — dO^ + sin^ 9 dtj)^. The lapse and shift vector are therefore: 



, r + 2M r + 2A/ 

The basic ADM variables read 

7„ = 1 + 2— ; = r' = 

r sm 6 

Krr = -^{r + M)a ; Kee = 2Ma = 



r 2Af 



sm^ ' 



and the CADM variables 



7r 



Iwr/ 3 • 2/>ni/3x 2Ma{r + 3M) 

-ln(r(r + 2Af)r'^sm^6l ' ) ; K ^ — —-, 

4 ' ^ r2 (r + 2M) 

r + 2M _ _ r _ 



r2 [{r + 2M) sin^ 6] ^'^ ' [{r + 2M) sin^ 6] ^'^ sin^ ( 

4M a{2r + 3Af) 
3 r4 [(r + 2Af) sin^ 6] ^'^ 
2M a(2r + 3Af) 



3 r(r + 2Af ) [(r + 2A//) sin^ O] ^'^ sin^ 6* 
4r3(r + 3Af)sin^6l 



3 r2(r + 2Af)5/3 
2 r(r + 2A/) cos 6* 
"3r2(r + 2Af)2/3sin2/3 6» 
f'^ 0. 

Note that some of the quantities are functions of 6*. In our spherically symmetric implementation of these equations 
we have explicitly expressed each variable as a function of r times the exact function of the angle 9. For instance we 
write 

^gg = hee{r)/sin^e. (4.3) 

Proceeding this way allows for the explicit appearance of 6 to drop out of the equations, providing at the end of the 
day, a truly ID system of equations corresponding to spherical symmetry. 



A. Comparison 

Extended tests were performed with both codes (under the same conditions) to understand the robustness of each 
formulation with excision. As has been observed in previous work CADM gives long term evolutions when the 

evolution of K is "frozen" ; ie the equation for K is not evolved or the value of K is fixed by the choice of a slicing 
that leaves K fixed (for instance, maximal slicing that fixes K = 0.). On the other hand, longer term evolutions have 
also been achieved with an area locking gauge in the ADM formulation p^. We then perform three basic tests: 
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• Fully free evolution: All equations corresponding to each system are integrated without imposing any further 
condition 

• "Locked" evolution: Conditions on some of the field variables are enforced (see below). 

• "Perturbed" evolution: Same as the "locked" case but considering linear perturbations of Schwarzschild space- 
time as initial data 

In all these tests, we study the dependence of the obtained solution under discretization size and location of the outer 
boundary. The inner boundary is placed at r — M and the outer boundary is varied (placed at r = nM) while 
keeping Ar = const. Outer boundary data are provided by 'blending' the numerical solution to the analytical 
one. This choice reduces gradients and second derivatives at the boundary allowing for a clean evolution without 
much reflections from the outer boundary. 



1. Fully free evolution 



In this case, all equations corresponding to systems (3.1,3.2) are evolved and the obtained solutions are compared 



We use the Hamiltonian and momentum constraint as monitors of the quality of the evolution. Our results can be 
summarized as follows. For the ADM formulation we observed that stable evolutions are obtained if n < 6 while for 
larger n > 6 the solution exhibits exponentially growing modes. It is worth emphasizing that the evolution is not 
unstable under the strict sense (i.e. the solution can be bounded from above by an exponential p^ .). However, the 
presence of this exponential mode will likely spoil any long term simulation. For the CADM system, irrespective of 
the value of n exponential modes are clearly present in the solutions. 

These results are illustrated in Fig.|^, which shows the L2 norm of the Hamiltonian and momentum constraints of 
the solutions obtained with both formulations when n = 4. Clearly, the solution obtained with the ADM is better 
than that obtained with the CADM. Figure ^, also displays the L2 norm of the Hamiltonian constraint for the case 
n = 9, although the solution obtained with the ADM formulation can be considered better than that from the CADM, 
both grow exponentially. 



A B 



1.0 ^ , ' , ' , ' , ' 1 6.0 

il 
'I 




0.0 50.0 100.0 150.0 200.0 250.0 0.0 20.0 40.0 60.0 80.0 

tiM] tm 

FIG. 2. The L2 norm of the Hamiltonian and momentum constraints for both formulations (discretization size Ar = M/10) 
for the case where the domain of integration is [M, 4M] (A) and [M, 9M] (B). In A, the evolution obtained with the ADM 
formulation does not show the presence of an exponentially growing mode like the one obtained with the CADM approach. 
However, for the larger domain, (B), solutions obtained with both formulations are exponentially growing. 



2. Locked evolution 

It has been observed in the literature that in the case where K is fixed in time very long term evolutions can 
be obtained with the CADM system. It is interesting to see what this condition implies, note that if K does not 
vary, then, (j) will remain unchanged and therefore the determinant of the three metric 7^^ will remain independent of 
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time. This, could be regarded as an evolution that "locks" the volume and bears some similarity with the so-called 
"area- locking" gauge jl^j2^. It is also worth pointing out that a similar strategy can be implemented in the ADM 
case as it has been shown in [|2j . In this work, by choosing a gauge that "locks" the evolution of gee the exponentially 
growing modes displayed by solutions in domains with 6 < n < 11 are removed. Figure ^ shows the L2 norm of the 
Hamiltonian and momentum constraints corresponding to solutions obtained with both formulations for the choice 
71 = 4 and n = 9. In both cases the simulations can be performed for unlimited time without observing exponential 
modes. Again, the solution obtained with the ADM formulation is slightly more accurate than that provided by the 
CADM formulation. 

It is important to observe that for the CADM case with K frozen, evolutions without exponentially growing modes 
are obtained with n as large as 16. For n > 16 long term evolutions (> lOOAf) display at late times a clear exponentially 
growing mode. 



0.20 



0.15 ^ 



0.10 ' 



0.0 



M„ 




100.0 



200.0 



0.06 t 



0.04 1 'J ^ 




0.0 



50.0 



100.0 



150.0 



200.0 



250.0 



FIG. 3. The L2 norm of the Hamiltonian and momentum constraints for both formulations (with Ar = M/10) for the case 
where the domain of integration is [M, 4M] (A). Neither formulation displays exponential modes in this domain and the ADM 
one yields more accurate results. In B, the L2 norm of the Hamiltonian and momentum constraints of the solution obtained 
with the CADM and an "area-locked" ADM evolution (in the [M, 9M] domain) is shown, a transient oscillatory behavior is 
present at earlier stages and the solutions then settle to a constant value. 



3. Perturbed evolution 



In this case, we test the evolutions under perturbations (using a locked evolution in the CADM case but not in the 
ADM one). The initial data corresponds to the analytic value of jrr (or ^rr) plus some arbitrary pulse of compact 
support. Of course, this data is unphysical but we use it to probe for stability of the implementations in a non-trivial 
scenario. The amplitude of the pulse is chosen such that it can be considered a linear perturbation of a Schwarzschild 
spacetime. The results obtained with both codes agree with those of the previous section. Figure |4| corresponds to 
the evolution of a pulse with compact support in [3M, 5M] being evolved in a computational domain of [M, 6M]. The 
L2 norm of the Hamiltonian constraint, after some initial transient behavior, settles into a stationary regime. 
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FIG. 4. The L2 norm of the Hamiltonian constraint for the perturbed evolutions (where K has been frozen in the CADM 
evolution while all fields are evolved in the ADM one). After some initial transient behavior both settle into an stationary 
solution. 



B. Locking evolutions 



As mentioned, in previous work by not evolving the equation for K very long term evolutions have been 
achieved with the CADM formulation. However, choosing to do so is unphysical in generic situations. One would like 
to have a prescription where a similar condition can be enforced without having to not evolve one or more equations. 
Here, one can use the gaug e free dom of the theory to demand that dtK = 0. This in turns implies a condition of the 
shift vector from equation ( 3.2c ), 



1 



(4.4) 



This condition is straightforward to implement in ID but is certainly more complicated in 3D. Additionally, there is 
a great deal of ambigui ty a s it is only one equation for three variables Therefore two supplementary conditions 
must be chosen so that (4^) can be used to 'freeze' t he ev olution of K. In our present implementaion we have simply 
chosen = (with A = {9, (f)) and obtain with (O) as 



1 



drK 



u ^3 



(4.5) 



A simple way to obtain is by a first order approximation of the rhs of Eq. [4.5| (ie. evaluating each term at the old 
level). By using this condition, instead of choosing not to evolve K, we were able to obtain evolutions not displaying 
exponential modes for times larger than 250A/ (with resolutions of Ar = M/IO and finer). Figure ^ illustrates what 
is obtained in a simulation with computational domain defined by [M, IIM] (with Ar = M/10). The values of the L2 
norms of the Hamiltonian, the function F — Ft=() and value oi K ~ Kt=o are shown as a function of time. Since is 
obtained only as a first order approximation, K and F are expected to vary during the evolution. As can be seen in 
Fig. H both grow linearly but stay fairly close to zero and the evolution proceeds without displaying an exponential 
growth. 
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FIG. 5. The L2 norms of the Hamiltonian constraint and the differences F ~ Ft = and K — Kt^o vs. time. The evolution 
proceeds without displaying exponential modes and the value of \K — Kt=o\ stays close to zero. 



C. Causal differencing and domain of dependencies 



As a last point, it is interesting to see how causal differencing is indeed providing a correct way to discretize the 
equations taking advantage of the causal properties of the spacetime. The fact that the null cones (or the causal domain 
of dependence) are tilted inside the horizon, allows for a stable numerical irnplementation where inner boundary data 
need not be provided if the inner boundary is inside the black hole (see Fig. g). This is possible because the numerical 
domain of dependence naturally contains the domain of dependence of the inner boundary point. This is of course, 
not true if the innermost point is outside the event horizon. In order to illustrate this fact we compare 2 cases (with 
both formulations) where the innermost point is placed inside or just outside the event horizon. As illustrated in 
Fig. 0, while the solution obtained with the inner boundary inside the black hole is stable, the other, as expected is 
unstable. 



P 

e 




EH 



FIG. 6. Domains of dependence of points inside and outside the event horizon (EH). Inside the EH, the past null cone of ; 
is tilted, therefore the evolution algorithm does not need the value of the fields at p on the old level. 
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FIG. 7. The L2 norm of the Hamihonian constraint for the cases where the inner boundary is inside (at r — 1.5M denoted 
with a sohd hne) and outside (at r = 2.2M denoted with a dashed hne) the event horizon (using the ADM formulation). 
Since the latter does not respect the CFL condition, the obtained solution is unstable. (The values shown correspond to a 
discretization of size Ar — M/10 and no qualitative difference is observed with finer resolutions.) 

V. CONCLUSIONS 

The results presented in this work show that excision techniques can be straightforwardly used in the CADM 
formulation directly from the structures developed for the ADM formulation. 

The ADM formulation is superior to the CADM both in accuracy and total time evolution when the evolution of 
K is not locked in CADM. When locking is implemented, then CADM is better than ADM as the solution obtained 
with the CADM formulation does not display exponential modes with the outer boundary placed as far as 16M. On 
the other hand, evolutions with the ADM formulation display exponential modes with the outer boundary placed 
at IIM and beyond. Additionally, for the case where outer boundaries are placed 'very' far, although exponentially 
growing solutions are present in solutions obtained with both formulations, ADM simulations crash at earlier times 
than CADM. 

It is worth remarking again that in both formulations, implementing a gauge that minimizes the changes in some 
of the fields (like gee in the ADM formulation or K in the CADM one) dramatically improves the evolutions in ID. 
In the 3D case, the use of 'area or circunference locking'^ is indeed more complicated than locking K, simply by the 
fact that in the former case one is trying to control a tensor component while in the latter a scalar. Thus, locking K 
is likely to have an easier and probably more general implementation than area locking (although in cases where the 
final black hole is close to a non spinning one, this implementation is rather straightforward). Controling the evolution 
of K demans a condition like that given by Eq (4.4), and two extra conditions on /3* will be required. An option that 



could mimic the ID implementation would be to foliate the 3D hypersurfaces with a sequence of 2-surfaces defined 
by 6 = const (with 6 the expansion of outgoing null rays). Once this foliation is obtained, the shift vector /3* could 
be decomposed as 



(5.1) 



with parallel (perpendicular) to the normal of the 2-surfaces. Thus, the two further conditions can be chosen 

such that I3\ = const. This will thus minimize changes in transversal directions and will resemble that we have used 
in the ID case. Of course, this is just one possible approach and further studies will be required to obtain a K fixing 
condition that leads to a practical implementation. 

In conclusion, implementing singularity excision techniques in the CADM formulation is straightforward. The 
usefulness of this implementation depends on implementing a gauge controlling the behavior of K . Assuming this can 



ie. controling the determinant of the angular part of the metric (area locking) or ggg (circunference locking) 
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be achieved, CADM appears to be capable of providing more robust simulations than ADM when the outer boundary 
is placed farther than IIM from the final black hole of mass M, if the outer boundary is closer, then the ADM 
formulation provides evolutions as stable as the CADM one but with better accuracy. 

As a last remark, it is important to stress that we have only implemented the causal differencing algorithm described 
in |2^Jl^ since at present is the only one fully implemented in 3D. Other alternatives have been proposed ]r^ , [T8| , p5t ; 
due to the restriction to spherical symmetry it is likely that the application of these will yield similar results to those 
presented in this work. 
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